Modeling relaxation and jamming in granular media 
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We introduce a stochastic microscopic model to investigate the jamming and reorganization 
of grains induced by an object moving through a granular medium. The model reproduces the 
experimentally observed periodic sawtooth fluctuations in the jamming force and predicts the period 
and the power spectrum in terms of the controllable physical parameters. It also predicts that the 
avalanche sizes, defined as the number of displaced grains during a single advance of the object, 
follow a power-law, P{s) ~ s""^, where the exponent is independent of the physical parameters. 
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Granular materials are composed of many solid parti- 
cles that interact only through contact forces, displaying 
a variety of behavior that distinguishes them from other 
forms of matter such as liquids or solids. While granu- 
lar materials can flow like a liquid, unlike liquids, they 
reach a jammed state when stressed In the jammed 
state, which has analogies in a variety of physical sys- 
tems such as dense colloidal suspensions, traffic flows, 
and spin glasses [§J|], the local dynamics of grains is 
frustrated by close contacts between neighboring grains. 
Although jamming in granular materials has previously 
been discussed in the context of the gravitational stress 
induced by the weight of the grains, it can be induced by 
any compressive stress, such as the stress generated by 
an object traveling through a granular material. Indeed, 
there is experimental evidence that a solid object be- 
ing pulled slowly through a granular medium is resisted 
by local jamming, and can only advance with large-scale 
reorganizations of the grains. The jamming and reorga- 
nization phenomenon, which can be detected through the 
drag force acting on the object in granular medium, re- 
flects both the nature of stress propagation through the 
medium and the dynamics of the granular medium [p|-p^ . 
The drag force opposing the motion of the object origi- 
nates in the force needed to induce such reorganizations 
and exhibits strong fluctuations with a stick-slip charac- 
ter associated with the reorganization of the grains d). 
While these phenomena were documented extensively ex- 
perimentally, a theoretical microscopic study of the be- 
havior of jamming and relaxation of the drag force has 
never been attempted. 

In this paper we introduce a microscopic model for the 
motion of a vertical cylinder through a granular medium, 
describing the jamming and reorganization of grains in 
terms of compression and relaxation of elastic springs 
| p6[ with random thresholds. The model can reproduce 
many of the experimentally documented features of the 
drag force, such as the sawtooth shape of the tempo- 
ral behavior and the main features of the power spec- 
trum. It predicts the period of the sawtooth pattern in 
terms of physical parameters such as depth, diameter, 



and cylinder velocity in agreement with recent experi- 
mental results 1^,^,0. In the jammed state, it predicts 
an avalanche-like relaxation of the grains, the avalanche 
size following a power-law distribution. The model also 
predicts that the critical relaxations generating the sud- 
den drops are nucleated preferentially from the bottom 
part of the cylinder. Finally, we find that the tempo- 
ral pattern of the drag force can change depending on 
the elasticity of the grain medium: when grains are rigid 
the temporal pattern is periodic and the avalanches fol- 
low a power-law, however, when they are more elastic, 
the temporal pattern is random and the avalanche size 
decays exponentially. 

Microscopic model — The model was constructed to 
emulate the drag experienced by a vertical cylinder of 
diameter d inserted to a depth in a granular bed [^,|| . 
The grains move with constant speed v in the positive 
X direction, pushing the cylinder in the same direction. 
The motion of the cylinder is constrained by a fixed stop 
which is coupled to the cylinder through an external 
spring with spring constant K. The force on the fixed 
stop is equivalent to the drag force, F[t), on the cylinder 
as a function of time. We refer to the spring located be- 
tween the cylinder and fixed stop as the external spring. 
As the granular medium moves in the positive x direction, 
the grains in contact with the cylinder's surface push the 
cylinder with small forces whose sum is the total drag 
force. To model the heterogeneous nature of granular 
drag we regard the surface of the cylinder as a planar 
rectangle partitioned into d x H cells of unit size. 

Grains push the cylinder only if they are compressed, 
and we thus model each cell as a "grain-spring". The 
magnitude of the force /j ^ (i) exercised by the spring in 
cell {i — l..d, j — 1..H) is given by Hooke's 

law, fi,j = fij{0) + kijXij where /i,j(0) is the x com- 
ponent of the ambient force, kij is the grain spring con- 
stant, and Xij is the deviation of the spring's length from 
its uncompressed value. It is well known that the pres- 
sure in granular media increases linearly with depth and 
will increase the ambient force /ij(0) = foj The 
grain-spring constant should be interpreted as describing 



the elasticity of the force chains instead of the individ- 
ual grains. The force chains are expected to get stiffer 
with depth, since the participating grains are more com- 
pressed, allowing less room for configurational changes. 
Thus we expect that the spring constant kij will also 
increase linearly, kij = koj. 

If the grains are too compressed, they will fail by slip- 
ping relative to the cylinder' surface and each other, thus 
relaxing the local forces. As a result, the total force act- 
ing on the cylinder decreases, and the cylinder slips rel- 
ative to the grains. To model this microscopic failure 
we introduce a critical threshold gt^, which is a random 
variable uniformly distributed between [goj, gij] , where 
go and gi are constants. When the elastic force on a grain 
spring exceeds its critical threshold, (fij > gi,j), the 
spring is relaxed to its equilibrium position, fij — /^^ (O), 
and the threshold gij is newly updated by a new random 
number. 

As grains advance in the positive a;-direction by the 
distance vSt during time interval dt, they increase the to- 
tal force acting on the cylinder, compressing the external 
spring K as well. The balance between the cumulative 
action of the grain springs and the opposing force of the 
external spring allows the cylinder to move in the -I- a; 
direction by a distance £, 

e = hv5t/ikt+K), (1) 

where kt = kodH{H + l)/2, is the collective spring con- 
stant of the grain springs. The distance (1) was ob- 
tained by balancing the collective elastic forces of the 
grain springs and the external spring on each side of the 
cylinder, 

J2h,jSx,,j =Ke, (2) 

where Sxij — vSt — £. After obtaining £, the effective 
compression of grain springs can be determined from 
Sxij = Kv5t/{kt + K), leading to an increase in the 
grain spring force by Sfij = kojKvSt/ {kt+K). When the 
grain springs are much stiffer than the external springs 
[kt ^ K), corresponding to the case in which the ex- 
periment was performed, the increased grain spring force 
acting on the cylinder at (i, j) becomes 

5fi,o = kojKvdt/kt. (3) 

The situation suddenly changes if a grain slips, i.e. the 
force fij on a grain spring reaches its threshold gij. We 
reset the force to fi,j{0), and the threshold is updated 
[ p^ . After this update, the balance between the elastic 
forces on each side of the cylinder breaks down, because 
the total force acting from the grains on the cylinder has 
dropped by fij{t) — /ij(0). As a result the cylinder will 
move backward (in the negative x direction), pushed by 
the external spring, compressing further the remaining 
grain springs. The displacement of the cylinder can be 
calculated by using the balance equation (2), where the 



newly updated sites (i.e., those having fij = fij{0)) are 
excluded from the summation. 

There are two possible outcomes of this slip event. 
First, if this sudden compression of all grain springs will 
not cause any more springs to reach their thresholds, af- 
ter establishing a new equilibrium we continue the con- 
tinuous compression of all springs by the motion of the 
grains with velocity v. However, in some cases the dis- 
continuous increase of the force on the grain springs will 
cause some other springs to reach their thresholds. In this 
case the updating (replacing each broken spring with an 
uncompressed one and calculating the new equilibrium) 
is repeated until no further reorganizations occur. The 
time is then incremented, followed by the advance of all 
grain springs, leading to a repetition of the above pro- 
cesses through compression and new updating. The dy- 
namics of the model are similar to the random fuse model 
in one dimension ]2C| ] describing fracture of a fiber bundle 
in the sense that when one spring (bond) breaks down, 
the load from it is shared by others. There is a significant 
difference, however, in that within the current model the 
springs can be re-compressed, resulting in a stationary 
dynamics, while in the random fuse model a bond is per- 
manently disconnected once it has burned out. 

The stochastic model described above offers a micro- 
scopic description of the system investigated experimen- 
tally. Despite its simplicity, as we show next, it accounts 
for many key factors of the observed behavior, and of- 
fers insight and quantitative predictions that were not 
available experimentally |^,^. 

Sawtooth pattern — A characteristic feature of the drag 
force observed in the experiments is that the force on a 
cylinder, F{t), increases linearly, followed by a a sudden 
drop in F{t), corresponding to a collective failure and 
reorganization of the grains. As Fig. 2a-c shows, this 
sawtooth pattern is fully reproduced by the model. The 
linear increase corresponds to a continuous compression 
of both the grain springs and the external spring. At 
a certain point, however, a grain spring fails, which re- 
sults in a collective and subsequent failure of all other 
springs in the system, since they are compressed to near 
their thresholds. Thus the stick-slip motion observed in 
the experiments correspond to two regimes: in the linear 
regime we see a linear convergence to the critical state, 
where all the springs are more or less simultaneously com- 
pressed towards their critical threshold. The sudden drop 
corresponds to an avalanche like spreading of a failure 
as soon as the critical or fragile state has been reached. 
The advantage of the presented model is that it allows 
us to quantitatively characterize the resulting stick-slip 
process. 

Linear Regime — What does determine the slope of 
the F{t) signal in the linear regime? Eq. (4) predicts 
that the drag force, F{t) — J^i j fi,j{t)j increases linearly 
with time with the slope ~ i^T in the jammed state. 
This linear increase is in complete agreement with the 
experimental results (see Fig. 2 of ||]). Furthermore we 
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predict that the slope is independent of the experimen- 
tal details, but depends only on the spring constant of 
the external spring, which is again consistent with the 
experimental findings. 

Failure and depth dependence — When updating oc- 
curs over the entire system, the drag force drops sud- 
denly, because most grain springs are reset to their equi- 
librium positions, and fi,j{t) = fijiO)- Accordingly, we 
expect the drag force to exhibit a sawtooth pattern. This 
is supported by extensive numerical simulations whose 
results are summarized in Fig.2a-c. This allows us to de- 
termine the average value of the drag force, that has been 
investigated extensively both experimentally and numeri- 
cally (see Ref. [5]). Since the force at is independent 
of i and proportional to j, we find that the average drag 
force F over time is proportional to ^ dH^. This result 
is confirmed by numerical simulation as well (see Fig.2d) , 
and is in agreement with the experimental results [p[. 

Power spectrum — To better characterize the fluctua- 
tions in the system, we measured the power spectrum 
of the time dependent drag force as predicted by the 
model. In Q it has been found that the power spec- 
trum is characterized by a few prominent peaks and sub- 
harmonics, determining the period of the signal, followed 
by a power-law tail at high frequencies which decays as 
f~^. The numerically determined power spectrum has 
the same features (Fig 2e), exhibiting an ~ 1//^ behav- 
ior at high frequencies and peaks at low frequencies. The 
origin of the 1//^ behavior comes from fluctuations of 
critical relaxation time due to random thresholds. The 
position of the largest peak corresponds to the inverse of 
the period of the sawtooth pattern. Thus, estimating the 
peak position, we can determine (l/T), the period of the 
sawtooth signal (see the inset of Fig.2e). Furthermore, 
this period can be predicted analytically as follows: the 
critical updating occurs when Sfij is increased to the 
maximum value of threshold jgi. The time required 
to reach this critical force through jamming is 

T ktigi - 9o)/{koKv) ^ (g, - go)H^d/Kv, (4) 

which represents the period. The numerical simulation 
data for different depths, diameters, and elastic spring 
constants K (shown in the inset of Fig.2e) confirm Eq.(4), 
and are also in agreement with recent experimental re- 
sults 

Avalanches — A quantity that has not been measured 
experimentally, but can be determined in the present 
model, is the avalanche size distribution P{s). When a 
collective failure occurs, this will result in the simultane- 
ous failure of a certain number of grains springs (but not 
all) creating an avalanche of failures. The avalanche size, 
s, is defined as the number of springs newly updated in a 
single advance of the cylinder. We find that the avalanche 
size follows a power-law distribution, P(s) ~ with 
T w 2.4(1) (see Fig. 3a). Since power-law distribution of 
events are common only at the critical point of spatially 
extended systems indicating that the continuous com- 
pression of the grains brings the system close to a critical 



state. Furthermore we find that the exponent r is uni- 
versal, independent of physical parameters such as depth 
H, width d, spring constant K, grain spring constant fco, 
and threshold gi — gg. Note that while the avalanche 
sizes could not be measured experimentally, bulk imag- 
ing techniques such as MRI, allowing one to compare the 
grain positions before and after an avalanche, could offer 
a quantitative check of our predictions. 

The simulations indicate that the critical avalanches 
are nucleated preferentially near the bottom of the cylin- 
der where the stress is largest (see inset to Fig. 3a). 
When a spring at large H is relaxed suddenly, the load 
taken over by the other springs is commensurately large, 
and has a higher probablity of nucleating a large-scale 
reorganization of the grains. In contrast, breakdown of 
springs at small H is less likely to nucleate avalanches. 

Finally, the model predicts that the clastic properties 
of the granular media have a strong impact on the tem- 
poral characteristics of the drag force. When grains are 
sufficiently soft, kt ^ K, the drag force develops a ran- 
dom signal rather than sawtooth, (inset to Fig. 3b). Such 
random characteristics also occur when the depth H is 
small enough to satisfy the condition kt ^ K for a given 
grain spring constant fcp. Interestingly, in this case, the 
avalanche size distribution does not follow a power-law, 
but decays exponentially as shown in Fig. 3b. 

In summary, we have introduced a stochastic model 
that describes the jamming and reorganization of grains 
associated with dragging an object through a granular 
medium. The model reproduces the sawtooth pattern 
of temporal evolution of the drag force and the 1//^ 
high frequency tail of the power spectrum, and predicts 
a power-law avalanche size distribution. This excellent 
agreement with the experiments is surprising because the 
model offers a mean-field treatment of the force chains 
which are known to be the basic mechanism of the stress 
propagation through grains. Indeed, the force chains are 
only included implicitly through the depth dependence 
of the grain-spring elastic constant. Thus improvements 
based on a more detailed handling of the force chains 
could be envisioned, but our model offers a crucial start- 
ing point for a detailed understanding of motion through 
granular media, and it offers a basis for more realistic 
modeling efforts. 
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FIG. 1. Schematic illustration of the stochastic spring model. The shaded area indicates the granular medium, which moves 
with velocity v in the positive x direction. The motion of the cylinder that tries to move along with the grains is opposed by a 
fixed stop coupled to the cylinder through a spring with spring constant K. We model the grains opposing the movement as 
springs with spring constant koj. 
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FIG. 2. (a-c) Plot of the drag force F{t) as a function of time for d = 2, H = 1000 (a), d = 1, H = 1000 (b), and d = 1, 

H — 500 (c), where /o — 0.5, go ~ 0.5, gi — 0.7, K — 1, ko — 1, and v = 1 arc used, (d) Double logarithmic plot of F/d versus 
H for different diameters d — 1 and 5. The solid line with slope 2.0 is obtained by a least square fit. (e) Double logarithmic 
plot of the power spectrum (PS) versus frequency. Inset: Double logarithmic plot of TK/d versus depth H for different spring 
constants K = 1 and K = 5 and cylinder diameters d = 1 and d = 5 of the cylinder. The data are well collapsed on the dotted 
line with slope 1.98, obtained by a least square fit, predicted by Eq.(4). 
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FIG. 3. a) Logarithmic plot of the avalanche size distribution P{s) versus size s for different diameters d = 1 and 2, depths 
H = 1000 and 2000, elgistic spring constants K = 1 and 10, and grain spring constants fco = 1 and 10, where critical avalanches 
which contribute an isolated point are excluded in the accumulation. All data, averaged over 5000 configurations, are collapsed 
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The inset shows the number of avalanches nucleated at a certain depth averaged over a 25 point depth interval. 



b) Scmi-logarithmic plot of the avalanche size distribution versus size for fco = 10 ^ showing an exponential distribution. The 
data are averaged over 5000 runs. The inset shows a plot of the drag force as a function of time under the same condition used 
in Fig.2b, but ko = 10"^ 
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